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Abstract. In this paper parametric formulas are 
presented to evaluate the flux of atmospheric muons 
in the range of vertical depth between 1.5 to 5 km 
of water equivalent (km w.e.) and up to 85° for 
the zenith angle. We take into account their arrival 
in bundles with different muon multiplicities. The 
energy of muons inside bundles is then computed 
considering the muon distance from the bundle axis. 

This parameterisation relies on a full Monte Carlo 
simulation of primary Cosmic Ray (CR) interactions, 
shower propagation in the atmosphere and muon 
transport in deep water [1]. The primary CR flux 
and interaction models, in the range in which they 
can produce muons which may reach ^1.5 km w.e., 
suffer from large experimental uncertainties. We 
used a primary CR flux and an interaction model 
able to correctly reproduce the flux, the multiplicity 
distribution, the spatial distance between muons as 
measured by the underground MACRO experiment. 

Keywords: Underground muon flux and energy 
spectrum. Simulation of muon bundles. Neutrino 
telescopes. 

I. Introduction 

The search for neutrinos of astrophysical origin with a 
km 3 scale neutrino telescopes is one of the most impor- 
tant challenges for high energy physics and astrophysics. 
Atmospheric muons represent the most abundant signal 
in a neutrino telescope, and are useful for calibration pur- 
poses and to check its response to the passage of charged 
particles. On the other hand, atmospheric muons may 
represent a background source because they may mimic 
high energy neutrino interactions; in particular, multiple 
muons in bundles seem to be particularly dangerous. A 
full Monte Carlo (MC) simulation of atmospheric show- 
ers can accurately reproduce the main features of high 
energy muons reaching a neutrino telescope, at a price 
of a huge amount of CPU time. Parameterisations of the 
underground/water atmospheric muon flux are available 
in the literature Q, [3); they allow an easy comparison 
of their predictions with experimental results. However, 
all of them do not take into account the simultaneous 
arrival of muons in bundles. 
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Fig. 1: Multiplicity distribution of muons in bundles 
evaluated with (Q3 at the depth of 3 km w.e. and for 
5 different values of the zenith angle. The computation 
assumes here a muon energy threshold of 1 GeV. 



The first author |2| assumes as muon integral energy 
spectrum at sea level the empirical formula given in (4). 
The author obtains an expression for the muon flux at 
depth >1 km w.e. and energy from few GeV up to 10 
TeV through a computation of the propagation of the 
particles in sea water. 

The second calculation [3| starts from an analytical 
expression of the muon spectrum at sea level and 
of the muon propagation in rock or water. The result 
is an analytical expression that gives the integral muon 
flux in the range 1 GeV -10 TeV for depths in the range 
1.5-16 km w. e.. The integral flux depends from 5 free 
parameters, which must be tuned in the expression of 
the sea level spectrum. 

In both cases no distinction is made between "single" 
muons and "multiple" muons (or muon bundles). These 
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fluxes are all-inclusive and contain a sort of summation 
over the multiplicity distribution, which actually changes 
with depth and zenith angle. 

Underwater/ice telescopes Q, Q, (U have a coarser 
spatial resolution with respect to underground detectors 
like MACRO. In particular, they cannot measure the 
muon multiplicity of atmospheric muon bundles. In addi- 
tion, their trigger and reconstruction efficiencies strongly 
depend from the bundle multiplicity. The present pa- 
rameterisation (full details in UJ) reproduces the flux of 
atmospheric muons as a function of their multiplicity, 
vertical depth and zenith angle (see Fig. [TJ. The energy 
spectrum is given as a function of the muon distance 
from the axis shower. 

From the parametric formulas, an event generator 
(MUPAGE, [9 1) was derived. It is well suited to be 
interfaced with MC simulations of atmospheric muons 
of underwater/ice experiments. Details and some appli- 
cation of this generator are presented in iflOl . 

II. From full MC to parametric formulas 

The present parameterisation of the multiple muon 
flux and energy spectra relies on a full MC simulation of 
the primary CR flux, interactions and shower propaga- 
tion in the atmosphere. The latest version of the HEMAS 
code was used (HEMAS -DPM flTJ). The reliability of 
the code is restricted to secondary particles with energies 
above Et = 500 GeV. Muons below Et (at sea level) 
are not followed through water because their survival 
probability for depths > 1.5 km w.e. is very small. The 
so called prompt muons (from the decay of charmed 
mesons) are not included. 

HEMAS was preferred to other codes since it was 
deeply used and cross-checked with the results of the un- 
derground MACRO experiment. Also the input primary 
CR spectrum was a phenomenological model obtained 
by MACRO (12). Both CR composition and some of 
the hadronic interaction features in HEMAS were tuned 
lfT3ll by studying the experimental muon multiplicity 
distribution [14|, the underground muon pair distance 
distribution (the so-called decoherence distribution lfT51 ) 
and the average energy of single and double muons lfl6l . 
These measurements are particularly important because 
constrain some quantities extrapolated in the MC out- 
side the kinematics region measured at accelerators, as 
the charged multiplicity distributions and the transverse 
momentum of muon parent mesons at TeV energies. 

The muons in the atmospheric shower reaching the 
sea-level (flat surface) were propagated down to 5 km 
of water using MUSIC fTTl . a 3D muon propagation 
code which uses recent and accurate cross sections of 
the muon interactions with matter. 

All the muon kinematics information from h =1.5 
km w.e. down to h =5.0 km w.e. were stored in ROOT- 
format files. We consider the vertical depth h as one of 
the main parameters, as the zenith angle 8. The flux is 
computed up to 8 = 85°. From these files the parameters 
of the analytic description of the muon flux and of the 
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Fig. 2: Muon vertical intensity versus depth of water or 
ice. Comparisons are shown among experimental points 
(reference in [3|), the ANTARES line 1 measurement 
lfT8l . our parameterisation (full line) and the one of [3| 
(dashed line). Both computations considered only muons 
with energy > 20 GeV at a given depth. The ANTARES 
measurement with 5 lines is reported in lfl9l . 



energy spectra were obtained with a fitting procedure. 
The energy of the muons in bundles depends (at a fixed 
h and 8) on the bundle multiplicity rn and on the distance 
R of each muon from the bundle axis. 



The flux (units: m s 



-2„-l„„-l 



) of muon bundles with 



multiplicity m (see Fig. [TJ is obtained as a function of 
vertical depth h and zenith angle 8 as: 

K(h,8) 



<P(m;h,8) 



jv(h,e) 



(1) 



The flux of bundles of increasing multiplicity m de- 
creases with increasing vertical depth and zenith angle. 
The parameters K(h, 8) and v(h 1 8) in (HJ depend from 
9 constants, which are reported in Q). From £Q, the 
depth-intensity relation for vertical muons (8 = 0°) can 
be calculated, and it is presented in Fig. [2] 

The energy spectrum of muons is described by: 

— ^— = GE.e^-^E, + e(l - e"<«)r 
d(log w E^) 

(2) 

where 7 is the spectral index of the primary CRs and 
e = a I (3 is defined from the parameters of the muon 
energy loss formula —{ dE jx^ ) = In this work, 
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7 and e were used as free fit parameters and (3 was fixed. 
The constant G = G(j, e) represents a normalization 
factor, in order that the integral over the muon energy 
spectrum (|2]l from 1 GeV to 500 TeV is equal to 1. In 
the case on single muon events (bundles with m = 1) 
the 7 = ~f(h), e = e(h,8) parameters depend on 6 
constants. 

Fig. [3] shows the evaluated energy distribution of 
single muons at 6 = and at various depths (h= 3, 4 and 
5 km w.e.). It is evaluated by multiplying the normalized 
energy distribution (|2]i with the flux ([T]) of single muons 
$i = $(m = 1, h,6 = 0). For instance, at ft, = 3 km 
w.e., $i = 1.6 x 10 -4 (to 2 s sr) _1 . The total muon 
flux at the same h and 9 is $t = Y2 m m &( m i h>> &) 
2.2 x 10~ 4 (m 2 s sr)^ 1 . For comparison, the result 
of the computation of J2) is reported. In this case, 
the total muon flux at ft = 3 km w.e., 9 = is 
2.1xl0 -4 (to 2 s sr)~ 1 . In order to compare with the 
distribution of our single muons, we normalized in Fig. 
[3] the prediction of [2] to the same number of muons. 
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Fig. 3: Expected energy spectrum of single muons at 
9 = and three different depths. The points represent the 
predictions (normalized to the same number of events) 
extracted from Fig. 3 of J2) (see text). 

The situation is more complicated for multiple muons. 
Due to the muon production kinematics, the muon 
energy depends from their distance with respect to the 
axis bundle. The description of the muon lateral distance 
R from the axis bundles is thus the preliminary step to 
evaluate the muon energy distribution in a bundle. R (in 
the plane orthogonal to the axis) was extracted from a 
distribution of the form: 



R 



— =C 
dR (R + R Q y 



(3) 



The parameters R a = R (h,m,9) and a = a(h,m) 
depend from 9 constants. 

The energy spectrum of muons arriving in bundles has 
the same general form as for single muons @. In the 
case of multiple muons, the analytic description of the 
parameters 7 = "f(h, R, to) and e = e(h,R,9) depend 
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Fig. 4: Differential energy spectra of double muons 
(to = 2), at five different distances from the bundle axis: 
R = 3, 10, 30, 50 and 70 to. The lines were computed 
with (O using the 15 constants reported in Tab. 6 of 
HI, assuming 9 = and h = 3.5 km w.e. = 3500 hg 
cm~ 2 . The marker points superimposed for references 
to the line R =3 m were obtained with the full MC 
simulation. 504 series of points were used by the fitting 
procedure to obtain the 15 constants. 



from 15 constants. An example of the energy distribution 
of double muons (to = 2) is shown in Fig. |4] 

III. Discussion of the results 

The present parameterisation allows to evaluate not 
only the total muon flux, but also the total number 
of muon bundles in deep detectors. The systematic 
uncertainties on the predictions are related to the all CR 
particle spectrum, chemical composition and interaction 
model. It should be noted that many experiments, like 
in neutrino telescopes, cannot measure the muon bundle 
multiplicity. The relevant quantity in calculating the total 
muon flux is the number of nucleons impinging at the 
top of the atmosphere. For the number of muon bundles, 
it seems more relevant the number of nuclei at the 
top of the atmosphere. The relative contribution to the 
total number changes significantly from one primary CR 
model to the other. 

The knowledge of the all-particle spectrum is limited 
by the gap from direct and indirect observations l22l . 
The overlap will be crucial in establishing a reference 
point of the latter and for modelling the origin of CRs 
at energies of the knee and beyond. At the energies 
of interest for neutrino telescopes, the primary CR 
composition deduced from experimental data shows a 
strong dependence on the interaction models adopted 
(see i.e. ||231 ). Deviations from one parameterisation to 
another vary from 10% to 100% in absolute value and, 
particularly, the position of the knee differs from one 
model to the other. 
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The L3+C data up to 2 TeV show that the muon flux 
predicted using different interaction models can differ 
by up to 30%. At higher energies the situation is worse. 
An overall uncertainty on the total muon flux and on 
the number of muon bundles vary from 30% to more 
than 100% at different depths and zenith angles. The 
measurements of the deep-intensity relation constrain the 
total vertical flux of muons at different depths (Fig. [2]i 
by not better than 30%. 

The approach of this work was to choose a combina- 
tion of the primary CR flux, CR composition and interac- 
tion model which reproduces (at the level of ^30%) the 
MACRO data (depth: 2000 - 5000 hg cm" 2 , 6 < 60°). 

Fig-E]shows the comparisons of the zenith distribution 
evaluated at a fixed depth using ([TJ and some experimen- 
tal data. 
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Fig. 5: Muon flux as a function of the cosO as measured 
by ANTARES OH, AMANDA-II |20|. and Baikal OH 
at three different depths. The results of our computation 
for a muon energy threshold of 20 GeV at the three 
depths is superimposed as full lines. 

The comparisons demonstrate that the analytical solu- 
tion provides a good representation of the experimental 
data so far collected. This is true also in the case of 
the Baikal experiment, at a depth nominally outside 
the range of validity of our parameterization. One has 
to consider that in many cases the analyses are still 
approximate and not all the experiments have declared 
their systematic uncertainties. 

The parametric formulas can be used also for under- 
ground experiments, provided that the correct conversion 
to equivalent standard rock is made. The parameters of 
the muon energy losses -( ^^ ) = a + [3E^ depend 
from the medium: at small depths, the water is a more 
effective absorber (per g/cm 2 ) because the ionization 
term a is greater than that for rock. At larger depths 
(> 1 km) the (3 terms begin to dominate and made rock 
a more effective absorber. This effect is illustrated in 
Fig. |6] 

IV. Conclusions 

A new parameterisation of the atmospheric muon flux 
under a homogeneous medium of water/ice, considering 
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Fig. 6: Factor by which depths in water must be multi- 
plied to convert to standard rock. For example, 1 .08 km 
of water are equivalent to 1000 hg cm -2 of standard 
rock Il24l . 



the simultaneous arrival of muons in bundles, is pre- 
sented. It allows to evaluate the total muon flux and 
the total number of muon bundles in deep detectors. 
Flux, zenith angle distributions and average energies 
obtained with these formulas are in good agreement 
with those predictions by other parameterisations and 
with experimental measurements. This parameterisation 
is well suited as a fast generator to simulate the atmo- 
spheric muon flux in underwater/ice neutrino telescopes, 
particularly for the neutrino background simulation. 
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